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Section  1 


INTRODUCTION 


Although  operational  retrievals  of  temperature  profiles  have  been  routinely 
performed  in  the  past  15  years,  the  techniques  and  procedures  used  have  been 
under  constant  refinement.  This  is  because  the  retrieval  of  atmospheric 
temperature  profiles  is  a  complicated  problem  in  which  consistent  and  accurate 
performance  is  difficult  to  achieve.  Many  retrieval  methods  have  been  proposed 
in  the  past;  these  can  be  categorized  into  two  groups:  statistical  and  physical 
methods . 

Statistical  methods  determine  temperature  profiles  based  on  results  from 
regression  analyses  of  collocated  satellite  and  sounding  data.  While  statistical 
methods  are  computationally  fast,  the  accuracy  of  the  retrieval  temperatures 
varies  significantly.  Physical  methods  determine  temperature  profiles  by  solving 
an  ill-conditioned  and  ill-posed  Fredholm  integral  equation  of  the  first  kind. 
Because  of  the  unstable  characteristics  of  the  solution  process,  conventional 
physical  methods  prescribe  first-guess  profiles  (for  reference)  and  specify 
external  constraints  to  control  the  behavior  of  the  solution.  For  example,  the 
minimum  variance  scheme  (Fleming  et  al.,  1986;  Eyre,  1989)  utilizes  numerical 
weather  prediction  models  or  statistical  library  search  procedures  to  generate 
first- guesses  and  the  use  of  elements  of  error  covariance  matrices  to  control  the 
smoothness  of  the  resulting  temperature  profiles. 

Current  efforts  focus  on  the  search  for  better  first-guess  profiles. 
Recently,  McMillin  (1991)  employed  a  classification  method  to  obtain  reasonable 
first-guess  temperature  profiles,  but  Improvement  of  the  retrieval  accuracies  was 
nominal  (-0.2-0. 3  K,  rms) .  Other  efforts  have  also  produced  inconclusive 
results.  In  the  BUAN  (Baseline  Upper  Atmospheric  Network)  experiment,  it  was 
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theorized  that  better  collocation  would  lead  to  better  retrievals .  For  this 
reason,  the  size  of  temporal  and  spatial  collocation  windows  were  reduced  from 
±3  hour  and  ±300  km  to  ±1  hour  and  ±150  km,  but  no  significant  retrieval 
improvements  were  achieved  (Bloom  and  Nappi,  1990). 

It  has  been  suggested  that  increasing  the  number  of  sensing  channels  may 
yield  additional  Information  content  for  temperature  retrievals  and  improve 
retrieval  accuracy.  Along  this  line,  simulation  studies  based  on  the  Advanced 
Moisture  and  Temperature  Sounders  (AMTS)  showed  that  better  retrievals  can  be 
achieved  by  using  high  resolution  channels  (Phillips  et  al. ,  1988).  Smith  et  al 
(1990)  demonstrated  that  root -mean -square  temperature  errors  can  be  reduced  to 
1-1.5  K  by  using  hundreds  of  channels  on  the  High  Resolution  Interferometer 
Sounder  (HIS)  while  onboard  aircraft,  but  only  limited  cases  for  a  short  time 
interval  and  over  a  small  geographical  locality  have  been  tested. 

In  an  attempt  to  circumvent  conventional  techniques,  King  (1985)  proposed 
a  novel  linear  transform  approach  to  solve  the  temperature  retrieval  problem. 
This  approach  is  referred  to  as  the  differential  inversion  method  (DI).  The 
essence  of  the  DI  is  that,  based  on  the  convolution  theoiry,  the  Laplace  transform 
of  Planck  intensity  profile  is  linearly  related  to  the  Laplace  transform  of 
radiance  profile.  By  applying  the  inverse  transform,  the  local  Planck  intensity 
can  be  expressed  as  a  linear  combination  of  the  derivatives  of  radiances  in  the 
logarithmic  peak  pressure  coordinate.  Instead  of  using  external  constraints  to 
"regularize"  the  retrieved  temperature  profile,  it  is  assumed  that  the  radiances 
can  be  represented  by  smooth  polynomial  functions,  so  that  the  derivatives  of 
radiances  can  be  computed.  Thus  the  DI  offers  an  alternative  approach  for 
solving  the  temperature  retrieval  problem  in  that  prescription  of  the  first- guess 
profile  and  the  associated  error  covariance  matrices  are  not  required. 
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Many  theoretical  aspects  of  the  DI  have  been  explored  (Llou  and  Ou,  1988; 
Liou  et  al.,  1989;  Ou  and  Llou,  1989).  It  has  been  shovm  that  the  DI  is 
efficient,  stable,  and  particularly  practical  for  narrow  weighting  functions. 
The  purpose  of  this  work  is  to  further  examine  the  applicability  of  DI  to 
temperature  retrievals.  Fundamentals  of  the  DI  are  presented  in  Section  2. 

We  have  acquired  the  BUAN  data  set  to  test  whether  the  DI  can  be  applied 
to  the  retrieval  of  atmospheric  temperature  profiles  using  satellite  data.  This 
data  set  contains  a  large  number  of  collocated  satellite  and  radiosonde 
observations  which  were  collected  between  January  15  and  July  15,  1988.  It  also 
Includes  results  of  temperature  retrievals  using  the  operational  method  for  both 
clear  and  partly  cloudy  conditions.  We  have  used  the  data  set  for  the  following 
purposes:  (1)  to  perform  simulation  studies  for  the  development  of  a  multi¬ 
channel  cloud-removal  method,  (2)  to  develop  an  adjustment  scheme  for  measured 
radiances,  and  (3)  to  compare  retrieval  results  using  the  DI  operational  methods. 
Details  of  the  BUAN  data  set  and  the  transmittance  model  for  radiative  transfer 
calculations  are  given  in  Section  3. 

To  apply  the  DI  to  retrieving  temperature  profiles  in  partly  cloudy 
atmospheres,  clear  radiances  must  be  reconstructed  from  cloud  contaminated 
radiances.  Many  cloud-removal  methods  have  been  proposed  in  conjunction  with 
various  retrieval  schemes.  Eyre  and  Watts  (1987)  have  reviewed  cloud-removal 
methods  which  have  been  developed  up  to  1987.  Most  of  these  methods  follow  the 
principles  of  the  N*-method  proposed  by  Smith  (1968).  The  quantity  N*  is  the 
ratio  of  effective  cloud  covers  for  adjacent  pixels  and  is  the  basic  parameter 
used  to  reconstruct  clear  radiances.  Different  approaches  have  also  been  used 
to  determine  N*,  which  are  basically  statistical  in  nature.  For  example,  Smith 
and  Woolf  (1976)  combined  HIRS  and  SCAMS  data  from  Nimbus  6  to  determine  N*  and 
clear  radiances  simultaneously  using  the  eigenvectors  of  the  covariance  matrix 
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of  clear  radiances.  McMillin  and  Dean  (1982)  enployed  radiances  of  both  HIRS  4.3 
fim  channels  and  MSU  (microwave)  channels  to  compute  N*  using  an  interchannel 
regression  technique.  We  have  developed  a  new  approach  that  uses  the  HIRS  Chs . 
6,  7,  and  8  radiances  to  simultaneously  determine  surface  temperature,  clear 
radiances  and  N*  based  on  results  from  radiative  transfer  simulations.  The 
method  solves  two  sets  of  equations  by  iteration.  One  set  is  the  multispectral 
correlations  of  the  simulated  clear  radiances  while  the  other  is  the  linear 
relationships  between  the  partly  cloudy  radiances  over  adjacent  pixels.  Details 
of  the  cloud  removal  method  are  provided  in  Section  4. 

While  the  DI  is  a  linear  inversion  method,  the  problem  involving 
atmospheric  temperature  retrievals  using  satellite  data  is  highly  nonlinear.  To 
apply  the  DI  for  temperature  retrievals  using  either  original  or  cloud-removed 
clear  radiances,  these  radiances  must  be  adjusted  so  that  major  nonlinear 
components  are  removed.  For  the  minimum-var lance  method,  adjustment  processes 
are  implicitly  included  in  the  specification  of  the  first  guesses,  the  error 
covariance  matrices  and  the  kernel  matrix  characterizing  forward  radiative 
transfer.  In  the  present  work,  adjustment  takes  the  form  of  correlation  between 
simulated  and  linear  radiances.  Discussion  of  the  adjustment  of  radiances  is 
given  in  Section  5.  Retrieval  results  and  discussions  are  presented  in  Section 
6.  Finally,  conclusions  are  presented  in  the  final  Section. 
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Section  2 


FUNDAMENTALS  OF  THE 
DIFFERENTIAL  INVERSION  METHOD 


2.1  Forward  Radiative  Transfer  Equation 

The  upwelling  radiance  at  the  top  of  the  atmosphere,  R„,  for  a  given 
spectral  band  may  be  derived  from  the  solution  of  the  basic  radiative  transfer 
equation.  Assuming  a  plane -parallel  atmosphere  under  local  thermodynamic 
equilibrium,  we  have 

R.  -  B,(T,)  r^(p.)  B,(p)  dp,  (2.1) 

where  is  the  Planck  intensity  at  central  wave  number  i/,  is  the 

transmittance,  p,  is  the  surface  pressure,  and  T,  is  the  surface  temperature. 
The  first  and  second  terms  represent  surface  and  atmospheric  emission 
contributions,  respectively. 

The  slab  below  the  surface  may  be  viewed  as  a  semi -infinite  isothermal 
emitter  whose  temperature  is  T, ,  and  whose  pressure  Increases  from  p,  downward 
to  «.  For  p  -»  ®,  r„(«)  is  taken  to  be  0.  Thus  r„(p,)  can  be  expressed  by  the 
following  mathematical  identity: 

r.(p.)  =  r,(p.)  -  r,(«)  -  J/-  dp.  (2.2) 

Multiplying  both  sides  of  Eq,  (2.2)  by  B„(T,)  ,  we  obtain 

B.(T.)  ^.(p,)  -  J;-B,(p,)  i.5^  dp.  (2.3) 

Combining  Eqs.  (2.1)  and  (2.3)  leads  to  the  following  forward  radiative  transfer 
equation  in  continuous  form; 
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•  r  mp)  dp,  (2.d) 

where  Bj,(p)  -  B^(T,)  for  p,  <  p  <  *. 

We  next  define  the  weighting  function  as  follows: 

tf,(p).-pf^.  (2.5) 

The  weighting  function  controls  the  relative  contribution  of  atmospheric  emission 
from  various  levels.  Because  r^(p)  decreases  as  p  increases,  aT„(p)/3p  <  0,  but 
*’^v(p)  ^  0  for  0  <  p  <  «.  Quantitatively,  3r„{p)/3p  depends  on  the  vertical 
distribution  of  absorbing  gases  and  temperature.  Near  the  top  of  the  atmosphere, 
rj,(p)  decreases  slowly  downward  due  to  the  low  density  of  absorbing  gases.  As 
the  density  of  absorbing  gases  increases  with  p,  the  rate  of  decrease  in  r„(p) 
increases.  Finally,  as  p  approaches  p, ,  T^,(p)  approaches  0,  and  ^T^,(p)/^p  also 
approaches  0.  Thus,  a  maximum  value  of  |aT„(p)/3pj  or  must  exist  between  the 
top  of  the  atmosphere  and  the  surface,  signifying  the  maximum  rate  of  change  of 
T^{p) .  Let  the  pressure  xevel  corresponding  to  maximum  be  p,  this  value 
varies  with  the  concentration  of  absorbing  gases,  temperature,  and  wave  number 
u.  If  the  wave  number  of  a  channel  is  near  the  band  center,  Tj,(p)  decreases 
rapidly  downward  due  to  strong  absorption,  so  that  p  is  located  at  a  relatively 
high  altitude.  If  the  wave  number  is  near  the  band  wing,  the  atmosphere  will  be 
relatively  transparent  and  r^(p)  does  not  vanish  at  p  -  p, ,  so  that  p  will  be 
near  the  surface. 

2.2  Differential  Inversion  Method 

We  shall  approach  the  inverse  problem  by  a  transformation  of  variables. 
The  following  variables  are  Introduced:  ^  -  -£np  and  ^  -  -?np,  where  p  is  the 
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peak  pressure  of  the  weighting  function  W.  In  terms  of  these  new  variables,  Eq , 
(2.4)  can  be  rewritten  as  (King,  1985;  Liou  and  Ou,  1988) 

R(?)  -  J‘  B«)  W((  -  Od^.  (2.6) 

where  R(?)  -  R^(p)  .  B(^)  -  B^(p),  and  W(^  -  ?)  -  W(p,p).  Assume  that  R  is  some 
smooth  function  of  the  variable  for  which  higher-order  derivatives  exist. 
This  assumption  implies  that  for  each  value  of  ^ ,  there  is  a  corresponding 
weighting  function  W(^  -  ^).  Subsequently,  we  perform  a  bilateral  Laplace 
transform  (Widder,  1971)  on  both  sides  of  Eq .  (2.6),  and  use  the  convolution 
theory  to  obtain 

r(s)  -  b(s)w(-s) ,  (2.7) 

where  s  denotes  the  transform  variable,  and  r(s),  b(s)  and  w(-s)  represent  the 
bilateral  Laplace  transforms  of  R,  B,  and  W,  respectively.  Based  on  Eq.  (2.7), 
we  can  obtain  6(C)  as  follows: 

B(?)  -  L-i 

[w(-s) 

where  L~^  is  the  inverse  bilateral  Laplace  transform.  The  term  l/w(-s)  may  be 
expanded  into  a  MacLaurin  series  (Pearson,  1971)  so  that 


(2.8) 


■  Is  ■ 


(2.9) 


where  the  coefficient  is  related  to  the  nth  derivative  of  the  function  l/w(-s) 
evaluated  at  s  -  0  as  follows: 


■  1  r 


/n! 


(2.10) 


By  inserting  Eq.  (2.9)  into  Eq.  (2.8),  and  noting  that  A„'s  are  constants,  we 
have 
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(2.11) 


B(T)  •  E  ts"r{s)]. 

n*0 


We  may  use  expressions  for  the  Laplace  transform  of  the  nth  derivative  of 
the  fvmction  R(4)  to  determine  the  inverse  transform  of  s"r(s).  Using  this 
procedure,  B(^)  can  be  expressed  in  terms  of  a  linear  sum  of  radiance  derivatives 
at  I  as  follows; 


B(^) 


V  d"R(0 

■*n  - — -  • 

AC 


(2.12) 


In  order  to  evaluate  higher  derivatives  on  the  right  side  of  Eq .  (2.12),  R(?) 
must  be  a  smooth  function  of  ^ .  In  principle,  the  solution  of  B(f)  would  be 
mathematically  exact  for  infinite  summations.  However,  for  narrow  weighting 
functions,  A„  converges  rapidly  to  zero,  so  that  the  summation  of  Eq.  (2.12)  may 
be  truncated  to  just  a  few  terms. 


2.3  Weighting  Functions  and  Inversion  Coefficients 

The  determination  of  B(f)  from  Eq.  (2.12)  requires  knowledge  of  both  and 
the  derivatives  of  R(f).  The  inversion  coefficient  A„  depends  on  the  weighting 
function  W.  The  bilateral  Laplace  transform  of  W  can  be  expressed  by 

w(-s)  -  J*p-*W(p)  df5/p,  (2.13) 

where  p  -  p/p.  In  principle,  if  the  atmospheric  temperature  and  humidity 
profiles  are  known,  it  is  then  possible  to  evaluate  w,  and  subsequently  l/w(-s) 
and  its  derivatives  numerically  through  accurate  line-by-line  spectral 
Integrations .  However,  because  atmospheric  parameters  are  not  known  before 
retrieval,  the  above  numerical  approach  is  not  practical,  and  we  prefer  to 
develop  an  analytic  form  that  can  adequately  approximate  the  weighting  functions. 
A  generalized  weighting  function  has  been  proposed  by  King  (1985)  in  the  form 
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W*<P)  -  «<*-!/«>  r-i(iA)p  exp(-pv«) . 


(2.U) 


where  r  is  the  Gasuaa  function  and  k  is  an  index  controlling  the  sharpness  of  the 
weighting  function.  Figure  1  shows  the  functional  form  of  for  x  -  0.5,  1,  and 
2.  Smaller  values  of  r  correspond  to  broader  weighting  functions  with  lower  peak 
values.  The  range  of  significant  contribution  of  the  weighting  function  is  from 
p  -  0.1  to  p  ••  10,  spanning  two  orders  of  magnitude. 

Using  the  functional  form  of  Eq.  (2.14),  we  can  show  that 

w(-s)  -  r[(l-s)/R)]/[r(l/R)R*/*]  .  (2.15) 

Evaluation  of  w(-s),  and  therefore  of  A„,  based  on  Eqs.  (2.13)  and  (2.10) 
Involves  digamma  and  polygamma  functions  (Abramowitz  and  Stegun,  1972). 
Expressions  for  as  functions  of  r  have  been  given  by  King  (1985)  for  n-1  to 
5.  Figure  2  shows  the  inversion  coefficients  X„  for  n-1  to  5  as  functions  of  r. 
For  K>1,  -1<A„<0.2,  and  A„  decreases  in  magnitude  as  r  increases.  In 
addition,  based  on  Eq.  (2.12),  any  error  inherent  in  the  derivatives  of  R(C) 
would  be  substantially  suppressed  in  the  retrieved  B(C),  because  A's  converge  to 
zero  for  larger  n.  For  r  <  1,  |Ai|,  lAg],  and  lAjj  are  still  less  than  1,  but 
IA2I  and  |A^|  Increase  sharply  with  decreasing  r.  Errors  in  radiances  will  then 
be  amplified  through  A2  and  A^  and  propagate  into  the  computation  of  B(^)  . 

2.4  Polynomial  Fitting  of  Radiance  Profiles 

Evaluation  of  higher  derivatives  of  R({)  Is  based  on  the  assumption  that 
radiances  can  be  represented  by  smooth  functions  in  terms  of  the  variable  We 
propose  to  express  R(^)  in  terms  of  an  (N-l)th  degree  polynomial  of  (  in  the  form 

Substituting  Eq.  (2.16)  into  Eq.  (2.12)  leads  to 
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d/d 


Fig.  1  The  generalized  weighting  function  as  a  function  of  p/p  for  the 

sharpness  index  «  -  0 . 5 ,  1 ,  and  2 . 
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Inversion  Coefficient 


Fig.  2 


Inversion  coefficients  Xj  (i  -  1.-5)  as  functions  of  the  index  k. 
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(2.17) 


_  *-1  ii-i  _ 

Equation  (2.17)  Is  also  an  (N-l)th  degree  polynonlal.  In  practice,  there  are 
only  a  finite  number  of  radiances  that  correspond  to  a  finite  number  of  The 
unknown  polynomial  coefficients  in  Eq.  (2.16)  can  be  obtained  by  least-square 
fitting  of  the  finite  number  of  R(fi)  to  a  polynomial.  Since  errors  in  the 
fitting  radiances  at  are  minimized,  the  first  term  of  Eq.  (2.12)  at  these 
would  be  more  reliable  than  at  other  f  values  between  neighboring  Note  that 
emission  from  regions  near  p^  makes  the  dominant  contribution  to  upwelling 
radiances,  and  it  would  be  physically  meaningful  to  present  results  in  terms  of 


the  temperatures  at  discrete  f values . 
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Section  3 


DATA  SET  USED  IN  TEMPERATURE  RETRIEVAL  STUDIES 


We  have  acquired  from  NESDIS  a  large  data  set  of  collocated  temperature 
profiles  and  satellite  data.  These  data  were  collected  during  the  observation 
period  (Jan  15  -  Jul  15,  1988)  of  the  Baseline  Upper  Air  Network  (BUAN) 
experiment.  Since  primary  goal  of  the  experiment  was  to  test  the  hypothesis  that 
the  accuracy  of  temperature  retrieval  can  be  improved  if  better  collocation  data 
are  available,  the  BUAN  data  set  represents  the  best  collocated  satellite  data 
set  so  far  available  for  temperature  retrieval  studies.  We  were  also  provided 
with  the  transmittance  model  currently  being  used  in  NESDIS  which  computes 
atmospheric  transmlttances  due  to  the  absorption  by  water  vapor  continuum  and 
rotational  bands,  CO2  15  /im  band,  ozone  9.6  fim  band  for  the  19  HIRS  channels. 

3 . 1  BUAN  Data  Set 

3.1.1  Satellite  Data  Set 

The  BUAN  data  set  contains  satellite  data  measured  by  the  radiometers  on 
board  NOAA-10.  The  data  of  NOAA-9  were  not  available  during  the  field  experiment 
due  to  an  Instrument  malfunction.  NOAA-10  had  local  equator  crossing  times  of 
approximately  0730  and  1930  LTC.  Except  in  polar  regions,  satellite  overpass 
generally  coincides  with  the  0000  and  1200  UTC  routine  synoptic  observations  at 
about  100  E  and  80  W  and  with  the  0600  and  1800  UTC  standard  observations  at 
about  30  E  and  170  W. 

There  are  two  levels  of  data  set.  The  level  1-T  (match  file  B)  contains 
original  brightness  temperature  data  for  the  19  HIRS,  four  MSU  and  three  SSU 
channels.  These  brightness  temperature  data  have  been  converted  from  counts 
recorded  by  the  radiometer  for  all  weather  conditions,  including  clear,  partly 
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cloudy  and  overcast  situations.  Each  record  contains  the  sounder  data  over  a 
"box"  which  consists  of  an  array  of  9x11  spots,  as  shown  in  Fig.  3.  Only  the 
data  for  the  9x7  spots  within  the  solid  fraae  of  Fig,  1  have  been  used  for 
retrieval  analysis.  One  or  more  radiosonde  stations  are  also  located  within  the 
domain  of  the  9x7  spots  ("+"  in  the  diagram).  The  BUAN  retrieval  procedures  have 
been  documented  by  Reale  et  al.  (1986).  A  brief  description  is  provided  in  the 
Appendix . 

The  other  level  of  data  (match  file  A)  contains  the  BUAN  retrieval 
products.  These  include  specifications  of  the  spot  for  which  retrievals  are 
performed,  the  associated  latitude, longitude,  date,  and  time,  the  "raw"  and 
cloud-removed  clear  radiances,  the  retrieved  temperature  and  humidity  profiles, 
the  average  value  of  N*  (used  in  the  cloud-removal  scheme),  surface  albedo,  skin 
temperature,  surface  elevation,  and  first-guess  temperature  and  humidity 
profiles.  It  is  noted  that  the  "raw"  clear  radiances  are  directly  adopted  from 
the  brightness  temperatures  for  an  identified  clear  spot,  and  that  the  "cloud- 
removed"  clear  radiances  are  generated  from  the  BUAN  N*-method. 

3.1.2  Collocated  Data  Set  of  Radiosonde  Observations 

BUAN  radiosonde  and  satellite  observations  were  collocated  and  saved  each 
day  during  the  period  of  data  collection.  The  general  procedures  were  to  collect 
pairs  of  radiosonde  and  satellite  observations  that  were  within  1  hour  and  150 
km  of  each  other.  [The  normal  operational  collocation  window  is  ±3  hours  and 
±300  km  (Kalnay  et  al.,  1990).]  The  radiosonde  data  contains  vertical  profiles 
of  pressure,  temperature,  dew-point  depression,  geopotential  height,  wind 
direction  and  wind  speed.  These  data  were  recorded  in  terms  of  15  mandatory  (all 
data)  and  25  significant  (temperature  and  pressure  only)  levels. 
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Data  spot 
Point  of  reti 
Radiosonde 


Display  of  the  array  of  11x9  spots  in  a  box.  Each  "x"  denotes  one  spot.  Cloud-clear  analyses  are 
performed  for  the  9x7  spots  within  the  solid  frame.  These  spots  are  further  divided  into  six  mini- 
boxes.  Only  one  retrieval  is  done  for  each  mini-box.  The  spot  over  which  a  retrieval  is  performed 
is  denoted  by  "o."  The  rawinsonde  station  is  denoted  by  "+."  The  distance  between  the  station  and 
the  retrieval  spot  is  less  than  150  km. 


3.2  Model  for  Computing  Transmlttances 

The  current  operational  transmittance  model  has  been  made  available  to  us 
through  NESDIS;  it  computes  the  vertical  profiles  of  transmlttances  involving 
TOVS  sounding  channels.  Spectral  bands  considered  Include  water  vapor  continuum 
and  rotational,  O3  -  9.3  na,  and  COj  -  15  /xm  bands.  For  the  water  vapor 
continuum  band,  effects  of  pressure  broadening  due  to  N2  and  water  vapor 
molecules  are  considered  based  on  physical  principles.  For  the  rest  of  the 
absorption  bands,  parameterizations  based  on  more  detailed  line-by-ltne 
computations  were  used  (McMillin  et  al.,  1979).  To  guarantee  adequate  match-up 
between  the  simulated  and  measured  radiances,  however,  an  adjustment  factor  7  was 
introduced.  This  factor  varies  with  wave  number  and  height.  A  single  7  is  used 
over  all  geographical  locations,  however.  Without  this  factor,  the  simulated 
radiances  will  differ  greatly  from  measurements  (Fleming  et  al.,  1986). 
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Section  4 


CLOUD-REMOVAL  METHOD 

We  have  developed  a  method  to  restore  clear  radiances  from  cloud- 
contaminated  brightness  temperatures  based  on  the  principle  of  the  N*-method 
(Smith  and  Woolf,  1976),  where  the  quantity  N*  is  the  ratio  of  effective  cloud 
covers  over  a  pair  of  adjacent  pixels.  The  current  method  differs  from  other  N‘- 
methods  in  that  N*  is  determined  using  the  radiances  involving  Chs.  6,  7  and  8 
based  on  radiative  transfer  simulations.  A  detailed  description  of  the  current 
cloud- removal  method  is  given  below. 

4.1  N*-Method 

For  a  pair  of  adjacent  pixels,  let  the  radiances  for  the  ith  channel  over 
clear  and  cloudy  areas  be  Roj  and  Rgi,  respectively.  Then  the  radiances  for 
pixels  1  and  2  are  given  by 

Ru  •  (1  -  fii)  Roil  ’?lReil»  (4.1a) 

Ri2  “  *  ^2^cl2’  (4.1b) 
where  and  are  the  effective  cloud  covers  for  the  respective  pixels.  These 
effective  covers  are  defined  as  the  product  of  cloud  cover  and  emisslvity.  We 
have  assumed  that  over  each  of  the  adjacent  pixels,  only  a  single  layer  of  cloud 
is  present,  and  that  clouds  are  located  at  the  same  altitude  with  the  same  cloud 
temperature.  In  addition,  clouds  are  assumed  to  have  the  same  radiative 
properties.  This  assumption  is  generally  valid  for  water  clouds  which  are 
thermally  black.  However,  it  cannot  be  applied  to  semi-transparent  cirrus 
clouds,  which  are  a  subject  currently  under  investigation.  The  above  assumptions 
lead  to  the  approximations  Ron  -  R^n  -  Roi  and  Ron  -  Rdz  -  Rd-  Consequently, 
we  can  eliminate  Rgi  from  Eqs.  (4.1a)  and  (4.1b)  and  obtain 
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N*  -  »r,/f?2 


(^.2) 


,  Roi  “  Rii 
^1  "  ^12 

If  we  know  N*  In  advance,  we  may  determine  from  Eq.  (4.2)  as  follows: 

Roi  -  (Rii  -Ri2N*)/(1  -N*).  (^-3) 

The  N*-method  is  subject  to  the  following  restrictions: 

1)  The  assumption  that  adjacent  clouds  are  of  the  same  height,  temperature, 
and  radiative  properties  may  not  be  realistic.  Clouds  may  occur  at 
different  heights  with  different  temperature  and  radiative  properties.  In 
addition,  multiple -layer  clouds  are  often  present.  If  the  cloud 
properties  for  the  adjacent  pixels  differ  drastically,  determined  from 
Eq.  (4.3)  will  be  subject  to  large  error,  and  therefore,  in  operational 
retrieval  methods,  efforts  have  been  made  to  search  for  the  pair  of 
adjacent  pixels  that  meet  the  assumption  most  closely  (McMlllin,  1978; 
McMillin  and  Dean,  1982) . 

2)  To  estimate  N*  from  Eq.  (4.2),  we  must  know  the  clear  radiance  of  a  channel 

which  is  not  used  in  temperature  retrieval.  This  Implies  that  the  surface 
temperature  must  be  known.  Over  partly  cloudy  and  overcast  areas, 

observed  IR  and  visible  radiances  are  contaminated  by  clouds  and  cannot  be 
directly  used  to  determine  surface  temperature.  An  alternative  approach 
is  to  use  microwave  channel  radiance,  but  microwave  measurements  are 
usually  subject  to  errors  caused  by  poor  resolution,  contamination  by 
precipitating  clouds,  and  uncertainties  in  the  estimated  surface 
emlssivity.  The  smallest  pixel  of  MSU  on  board  the  NOAa  polar-orbitxng 
satellites  has  a  diameter  of  about  110  km,  which  is  more  than  six  times 
larger  than  the  diameter  of  HIRS  pixels.  In  addition,  since  the 
wavelength  of  microwave  channels  is  on  the  ordev  of  mm,  the  measured 
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microwave  radiation  is  subject  to  scattering/absorption  effects  caused  by 
precipitation  drops  with  sizes  larger  than  or  comparable  to  the 
wavelength.  Finally,  the  surface  emissivity  in  the  microwave  region 
varies  over  a  wide  range.  Over  the  oceans,  the  emissivity  ranges  from 
0.45  to  0.65,  depending  on  temperature  and  surface  wind  speed.  Over  land, 
the  emissivity  ranges  from  0.90  to  0.95,  depending  on  soil  moisture.  Over 
ice,  the  emissivity  is  usually  >  0.7  (Burroughs,  1991).  It  is  generally 
recognized  that  applying  the  microwave  technique  over  land  and  ice  is 
extremely  difficult. 

3)  The  N*-method  cannot  be  used  when  both  pixels  have  the  same  amount  of  cloud 
cover.  In  this  case,  N*-l  and  the  denominator  in  Eq.  (4.3)  becomes 
infinite. 

4)  The  N*-method  will  also  fail  when  an  inversion  is  present  in  the 
temperature  profile.  In  this  case,  Ru  -  Ri2  “  ^i- 

4.2  Multi -spectral  Cloud  Removal  Scheme 

We  have  developed  a  multi -spectral  approach  using  the  HIRS  Chs .  6,  7  and 
8  radiances  to  determine  the  surface  temperature  and  clear  radiances  of  these 
channels  simultaneously.  First,  radiative  transfer  simulations  are  performed  to 
develop  correlations  between  clear  radiances  of  the  three  channels,  next  the 
relationships  are  examined  between  partly-cloudy  radiances  of  the  same  channels, 
and  finally  the  algebraic  equations  are  derived  and  solved  based  on  the 
simulation  analyses. 

4.2.1  Multi -spectral  Correlation  of  Clear  Radiances 

To  economize  the  computational  effort  involving  the  calculation  of  cl._ar 
radiances,  the  correlations  of  clear  radiances  have  been  developed  for  channels 
having  peaks  of  the  weighting  function  close  to  the  surface .  These  correlations 
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are  largely  independent  of  temperature  and  humidity  profiles.  The  BUAN  sounding 
temperature  and  humidity  profiles  have  been  used  to  perform  simulation  studies. 
HIRS  Chs.  6,  7,  and  8  are  selected  because  their  weighting  functions  peak  near 
the  surface  and  they  are  sensitive  to  surface  temperature  and  the  presence  of 
clouds.  The  correlations  are  developed  in  terms  of  ARj  -  R,i  -  R^i ,  where  Rqi  and 
R,i  are  the  clear  radiance  and  surface  emitted  radiance  for  Ch.  i. 

Figures  4a  and  4b  show  the  correlations  of  AR7  vs.  ARj  and  ARg  vs.  AR7, 
respectively.  The  quantities  R07  and  R^j  are  obtained  based  on  radiative  transfer 
simulations  using  a  large  number  of  temperature  and  humidity  profiles  contained 
in  the  BUAN  data  set.  Weighting  functions  are  computed  using  the  NESDIS 
transmittance  model.  The  correlation  between  Ch.  6  and  Ch.  7  radiance  difference 
is  extremely  good  because  the  surface  emission  contributions  have  been  removed. 
The  ARi  values  are  dependent  on  the  humidity.  Small  values  would  correspond  to 
polar  (dry)  conditions,  while  large  values  would  correspond  to  tropical  (moist) 
conditions.  For  AR^  >  40  mW/mVsr/cm"^,  the  computed  points  exhibit  larger  spread 
because  of  increasing  humidity  effects.  The  spread,  however,  is  within  only  3 
mW/mVsr/cm"^ .  The  solid  curve  in  Fig.  4  a  is  the  quadrature  fit  in  the  form 

4R, 

n»o 

where  d^  -  0.5763;  dj  -  0.3736;  d2  -  0.0047.  The  rms  difference  between 
simulated  and  fitted  values  of  AR7  is  about  1  mW/m^/sr/cm“^ .  Thus,  Eq.  (4.4)  is 
a  reasonable  representation  of  the  correlation  between  Ch.  6  and  Ch.  7  radiance 
differences.  On  the  other  hand,  the  correlation  between  Chs.  7  and  8  radiance 
differences,  as  shown  in  Fig.  4b  in  log-log  plot,  is  less  satisfactory.  For  AR7 
>  20  mW/mVsr/cm~^ ,  its  scatter  decreases  between  0.1  and  20  mW/mVsr/cm"^  as  AR7 
increases.  For  AR7  <  20  mW/mVsr/cm"^,  although  the  computed  data  scatter  widely, 
the  ARg  values  <  2  mW/mVsr/cm~^ .  Thus  it  is  reasonable  to  approximate 
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Continued. 


R,0  using  Rpg.  A  linear  regression  analysis  has  been  performed,  and  the  AR, 
values  are  correlated  with  AR7  as  follows: 

<  nARg  -  ♦  CjInAR;,  (4.5) 

where  e^  -  -14.33,  and  e^  -  4.35. 


4.2.2  Relationship  Between  Partly-Cloudy  Radiances  of  Selected  Channels 

The  N*-method  is  based  on  the  assumption  that  the  partly -cloudy  radiances 
for  a  particular  pixel  is  a  linear  function  of  the  cloud  cover  fj.  Given  that  N* 
is  constant  for  all  channels,  we  may  derive  from  Eq.  (4.2)  the  following 
equation: 


(4.6) 


In  the  two-dimensional  space  of  R^  and  Rj,  Gij  is  the  slope  of  the  straight  line 
connecting  the  points  (Rji,  Rii)  •  (Rj2i  R12) .  and  (Roj,  R^i).  Further  mathematical 
manipulations  of  Eq.  (4.6)  lead  to 


(4.7) 


The  unknown  R^^  and  Rqj  are  eliminated,  so  that  can  be  estimated  from  known 
values  of  R^i,  Ri2,  Rji  and  R_j2-  Once  is  obtained,  a  linear  relationship 
between  R^^  and  R^.  can  be  established  for  either  of  the  adjacent  pixels  based  on 
Eq.  (4.6).  This  linear  relationship  can  be  used  in  conjunction  with  the 


correlation  of  clear  radiances  to  solve  for  clear  radiances.  However,  Eqs.  (4.6) 
and  (4.7)  are  valid  only  when  the  assumptions  made  in  the  N*-method  are  valid. 


i.e.  ,  the  cloud  properties  are  the  same  over  adjacent  pixels.  Thus,  to  apply  the 


multi -spectral  method,  we  must  select  the  most  appropriate  value  of  Gij  among  a 
sample  of  adjacent  pairs.  This  is  described  in  the  next  subsection. 


23 


4.2.3  Estlnate  of  Gj^j 

Figure  5  illustrates  the  principle  Including  the  optimal  estimate  for  G^j . 
Shown  in  the  diagram  Is  a  two-dimensional  display  of  Chs.  6  and  7  radiances.  The 
inset  in  the  lower -right  portion  of  the  diagram  is  a  3x3  array  of  numbered 
pixels.  The  3x3  array  represents  one  basic  working  unit  of  the  cloud  removal 
scheme.  El^t  pairs  of  adjacent  pixels  are  formed  by  pairing  the  center  pixel 
(#9)  with  eight  surrounding  pixels.  The  nine  radiance  points  corresponding  to 
the  pixels  are  plotted  along  with  the  Planck  Intensity  curve.  These  radiances 
are  assigned  based  on  physical  principles.  They  are  smaller  than  the 
corresponding  clear  radiances  (Ro?.  Roe)-  Because  Ch.  7  is  more  transparent  than 
Ch.  6,  Ro7  >  Rq5,  all  partly-cloudy  Ch.  7  radiances  are  also  larger  than  Ch.  6 
radiances.  It  is  also  noted  that  since  the  wave  numbers  of  Chs.  6  and  7  are  very 
close,  the  Planck  intensity  curve  almost  coincides  with  the  curve  R^  -  R;.  The 
partly-cloudy  radiance  points  do  not  fall  on  a  single  straight  line.  The 
distribution  of  radiances  thus  signifies  a  tjrpical  situation  of  different  cloud 
properties  for  different  pixels.  This  situation  frequently  occurs  in  the  real 
data.  Ve  note  that  the  point  for  clear  radiances  is  denoted  as  "0",  which  also 
falls  on  the  curve  denoting  the  correlation  between  clear  radiances  (Section 
4.2.1).  The  clear  radiance  curve  intersects  the  Planck  intensity  curve  at  the 
point  (R,6,  R,?) .  which  corresponds  to  the  surface  emission.  Among  the  partly- 
cloudy  radiance  points,  point  #3  and  0  show  the  largest  difference.  The  slope 
of  93  (the  straight  line  connecting  i^9  and  #3)  is  close  to  that  of  90.  Using  the 
slope  of  93  as  G7(,  we  obtain  a  solution  for  the  clear  radiance  at  0'  (intersect 
between  93  and  the  clear  radiance  curve).  Evidently,  0'  is  different  from  0. 
However,  this  solution  is  more  accurate  than  the  solution  obtained  by  using  the 
slope  of  other  straight  lines  that  connect  point  #9  with  other  points.  Thus,  in 
practice,  we  select  the  pair  that  has  the  largest  difference  in  Ch.  7  radiances, 
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and  compute  its  associated  slopes  G74  and  Ggj.  This  strategy  does  not  guarantee 
the  accuracy  of  the  derived  clear  radiances,  however.  In  particular,  if  all  the 
radiance  points  are  close  to  each  other,  the  slope  determined  according  to 
this  strategy  may  still  have  a  large  error.  Thus,  additional  restrictions  are 
required  to  assume  that  the  determined  value  is  reasonable. 

We  have  performed  another  set  of  simulation  studies.  First,  the  cloud 
cover  and  cloud- top  pressure  is  prescribed  randomly  for  each  pair  of  (R„j,  R,,;)  • 
Next,  the  clear  radiances  were  computed  for  the  clear  radiance  correlation 
(Section  4.2.1).  The  cloudy  radiances  are  then  computed  based  on  the  following 
equation: 


Rci  “  B„(T,)r/p,)  +  JJ’B^(p)£^^dp.  1-6.  7.  8. 


(4.8) 


where  p,.  and  are  cloud- top  pressure  and  temperature,  respectively.  Partly- 
cloudy  radiances  are  obtained  from  Eq.  (4.1a).  Finally,  G74  is  obtained  from  Eq. 
(4.6)  using  the  simulated  clear  and  partly-cloudy  radiances.  (The  simulated 
values  of  G;^  and  G87  are  between  1.0  and  3.0.  Also,  G^j  increases  with  p^,) 
Based  on  this  analysis,  any  G^j  determined  from  Eq.  (4.7)  with  a  value  larger 
than  3.0  or  less  than  1.0  is  rejected. 

Finally,  we  combine  Eqs,  (4.4),  (4.5),  and  (4.6)  to  solve  for  clear 
radiances  of  Chs.  6,  7,  and  8  and  surface  temperature  simultaneously.  A 
numerical  Iterative  scheme  has  been  developed  for  this  purpose. 
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Section  5 


APPLICATION  OF  DI  TO  TEMPERATURE  RETRIEVAL 

5.1  HIRS  Channel  Characteristics 

In  order  to  Investigate  the  potential  applicability  of  DI,  HIRS  channels 
(Smith  et  al.,  1979)  are  used  to  perform  the  temperature  retrieval.  HIRS 
consists  of  19  channels:  Chs.  1-7  are  In  the  C02-15  ftm  band,  Ch.  8  Is  In  the  11 
urn  window  band,  Ch.  9  is  in  the  O3-9.6  fm  band,  Chs.  10-12  are  In  the  water  vapor 
vibration- rotation  band,  Chs.  13-17  are  in  the  C02-4.3  pim  band,  and,  finally, 
Chs.  18-19  are  in  the  3.7  nm  window  band.  For  the  purpose  of  temperature 
retrieval,  the  seven  channels  in  the  C02*15  pm  band  and  the  five  channels  in  the 
CO2-4.3  pm  band  will  be  used.  It  is  noted  that  the  weighting  functions  of  Chs. 
1-3  and  Ch.  17  peak  in  the  stratosphere  (p  <  100  mb). 

5.2  Parameters  of  Weighting  Functions 

We  use  the  generalized  weighting  function  defined  in  Eq.  (2.15)  in 
conjunction  with  retrieval  analysis.  The  key  parameters  in  this  function  are  k 
and  p.  Because  the  weighting  function  depends  slightly  on  atmospheric 
temperature  and  humidity  profiles  which  are  not  known  before  retrieval,  accurate 
values  of  k  and  p  cannot  be  computed  in  advance.  For  this  reason,  we  use 
constant  values  of  k  and  p  in  the  retrieval.  To  prescribe  k  and  p,  the  BUAN 
temperature  and  humidity  profiles  and  the  NESDIS  transmittance  model  are  used  to 
compute  a  large  number  of  weighting  function  profiles  for  HIRS  Chs.  1-7  and  Chs. 
13-17.  We  then  use  the  numerical  method  previously  developed  to  fit  each 
weighting  function  to  the  generalized  form  to  determine  k  and  p  (Ou  and  Llou, 
1989).  Subsequently,  we  average  all  ic  and  p  values  for  each  channel  to  obtain 
mean  values.  Figure  6  shows  the  average  values  of  k  and  p  for  each  channel 
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Distribution  of  the  mean  sharpness  index_ic  and  the  peak  pressure  p.  The  vertical  and  horizontal 
bars  denote  standard  deviation  for  k  and  p,  respectively.  The  numerals  denote  the  channel  number. 


along  with  the  standard  deviations  of  these  values  in  terms  of  horizontal  and 
vertical  bars.  It  is  evident  that  the  #c  values  for  Chs.  1-5  and  14-17  are 
smaller  than  2  with  standard  deviations  less  than  0.5.  In  contrast,  these  values 
for  Chs.  6,  7,  and  13  are  all  greater  than  2  with  standard  deviations  larger  than 
1.5.  Since  Chs.  6,  7,  and  13  are  affected  by  the  effects  of  water  vapor 
absorption/emission,  variation  in  k  is  primarily  associated  with  variation  in  the 
humidity  profile  near  the  surface.  The  range  of  p  is  generally  within  70  mb  of 
the  mean  values.  The  tropospheric  channels  have  larger  standard  deviations  in 
p  than  the  stratospheric  channels.  This  indicates  that  effects  of  variations  in 
p  are  not  negligible  for  tropospheric  temperature  retrievals.  To  incorporate  the 
effects  of  variation  in  k  and  p,  an  iterative  approach  is  required  to 
simultaneously  retrieve  both  humidity  and  temperature  profiles. 

5.3  Optimal  Degree  of  Polynomial  Fitting  of  Radiances 

Based  on  our  previous  investigation,  the  selection  of  the  degree  of 
polynomial  fitting  is  a  critical  factor  for  the  successful  application  of  DI  to 
temperature  retrieval.  Liou  and  Ou  (1988)  have  tested  the  accuracy  of  the 
retrieved  temperature  profiles  by  using  different  degrees  of  pol3rnomials  in  the 
synthetic  retrieval  exercises.  They  found  that  by  utilizing  the  radiances  of  the 
first  seven  HIRS  channels,  the  fifth-degree  polynomial  fits  almost  exactly  all 
the  radiances  and  yields  the  most  accurate  retrievals.  Figure  7  Illustrates  the 
principles  behind  the  selection  of  the  optimal  degree  of  polynomial  fitting.  The 
radiances  of  the  first  seven  HIRS  channels  are  fitted  to  polynomials  of  different 
degrees  (third,  fifth,  and  sixth).  Mathematically,  the  accuracy  of  the  fitting 
is  directly  proportional  to  the  degree  of  the  polynomial  used  in  the  fitting. 
The  sixth-degree  fitting  is  the  exact  fit  for  the  seven  radiances.  An  Important 
assumption  in  DI,  however,  is  that  the  functional  form  representing  the  radiance 
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Fitting  of  the  simulated  radiances  for  HIRS  Chs.  1-7  by  3rd-.  5th-.  and  6th-degree  polynomials 


profile  Bust  be  smooth.  From  this  figure,  it  is  evident  that  the  sixth- degree 
polynomial  shows  fluctuation  between  points  while  the  fifth-degree  pol3nnomial  is 
much  smoother.  In  view  of  the  preceding,  the  selection  of  the  optimal  degree  of 
the  polynomial  is  based  on  a  combination  of  the  smoothness  of  the  pol3momial  and 
the  accuracy  of  the  fitting.  If  the  number  of  channels  increases,  it  is  expected 
that  the  optimal  degree  of  polynomial  fitting  will  be  higher.  We  have  chosen  the 
fifth-degree  polynomial  as  an  initial  attempt  to  investigate  the  applicability 
of  DI  to  temperature  retrievals.  Our  previous  research  has  demonstrated  that  the 
fifth-degree  polynomial  fittings  lead  to  acceptable  retrieval  accuracy.  Further 
investigations  are  required  for  developing  a  systematic  method  to  determine  the 
optimal  polynomial  degree . 

5.4  Adjustment  of  Radiances 

Previous  synthetic  retrieval  studies  show  that  applications  of  the  DI  using 
simulated  HIRS  radiances  without  adjustments  produce  large  errors  (-5  K;  Ou  and 
Liou,  1989).  These  errors  occur  because  the  problem  of  remote  sounding  of 
temperature  is  nonlinear.  The  main  sources  of  nonlinearity  are  (Houghton  et  al.  , 
1984) : 

•  dependence  of  the  transmittance  on  temperature 

•  dependence  of  the  Planck  intensity  on  wave  number  within  the 

spectrum  of  a  sensing  channel 

•  dependence  of  the  Planck  intensity  on  wave  nvunber  between  spectral 
bands 

•  dependence  of  the  transmittance  on  the  concentration  of  absorbing 
gases  and  aerosols 

•  model  constraints 

•  effects  of  clouds 
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The  first  of  these  sources  is  relatively  weak,  but  other  factors  may 
significantly  affect  retrieval  accuracy.  We  have  already  described  a  method  to 
remove  cloud  contamination  in  satellite  data  in  the  previous  section.  Clear 
radiances  are  obtained  either  by  using  values  of  an  identified  clear  pixel  or  by 
applying  the  cloud-removal  procedure  involving  a  partly-cloudy  pixel.  To  obtain 
satisfactory  retrieval  results,  these  clear  radiances  must  be  "linearized"  in  the 
context  of  Dl.  Ou  and  Llou  (1989)  have  proposed  an  adjustment  scheme  in  terms 
of  a  "scaling  factor"  which  encompasses  the  effects  of  variations  in  the 
sharpness  index,  surface  discontinuity  of  temperature,  and  channel  properties. 
The  present  adjustment  process  is  based  on  this  schen.-  .  Furthermore,  we  develop 
a  correlation  between  original  and  "linear"  radiances.  The  linear  radiances  are 
generated  from  the  weighting  function  profile  with  uniform  shape  (fixed  k)  and 
the  Planck  intensity  profile  at  a  reference  wave  number,  subject  to  model 
constraints  that  are  explained  in  the  next  paragraph.  All  the  preceding  major 
nonlinear  factors  may  be  removed  by  the  procedure.  The  correlation  is  not 
sensitive  to  variations  in  temperature  and  humidity  profiles  (Ou  and  Liou,  1989)  . 

An  important  model  constraint  of  DI  is  that  it  is  developed  for  an  infinite 
atmosphere.  The  real  atmosphere,  however,  is  semi-infinite  with  the  surface  as 
the  lower  boundary.  In  Section  2,  we  have  demonstrated  that  the  radiative 
transfer  equation  for  a  semi-infinite  atmosphere  denoted  in  Eq.  (2.1)  may  be 
transformed  to  an  equation  corresponding  to  an  infinite  atmosphere  [Eq.  (2.4)]. 
The  basic  assumption  is  that  the  slab  below  the  earth  surface  may  be  viewed  as 
a  semi- infinite  isothermal  emitter.  This  conversion  implies  that  the  pol3momial 
expression  for  the  Planck  intensity  profile  [Eq.  (2.18)]  should  cover  the  whole 
range  of  pressure  (p  -  0  to  p  ®)  .  Discontinuity  always  occurs  between  the  skin 
temperature  and  surface  air  temperature.  Polynomial  expressions  generally  cannot 
adequately  represetit  this  discontinuity  and  the  behavior  of  constant  temperature 
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below  the  surface.  Figure  8  Illustrates  this  point.  The  solid  line  is  the  HIRS 
Ch.  7  Planck  profile  based  on  the  US  Standard  temperature  profile.  Below  the 
surface,  temperature  is  assumed  to  be  constant.  Also  shown  is  the  Ch.  7 
weighting  function  based  on  the  generalized  form.  The  dashed  line  s  the  fifth- 
degree  polynomial  which  fits  the  full  range  of  temperature  profile  for  30  mb  < 
p  <  200  mb  (W  >  0.02).  Clearly,  this  fitting  is  not  satisfactory.  Not  only  is 
the  below-surface  constant  temperature  not  able  to  be  simulated  well,  but 
significant  deviations  also  occur  for  the  stratospheric  profile.  The  fifth- 
degree  polynomial  fitting  is  much  improved,  but  only  for  the  atmospheric  portion. 
The  part  of  the  temperature  profile  below  the  surface  is  not  importanr  in 
retrieval  practice.  For  this  reason,  we  may  confine  the  fitting  to  the  above¬ 
surface  (p  <  1000  mb)  temperature  profile.  The  linear  radiance  is  evaluated 
based  on  the  polynomial  expression  that  closely  fits  the  atmospheric  {above - 
surface)  Planck  intensity  profile. 

To  correlate  the  satellite  radiances  and  linear  radiances ,  we  use  the  BUAN 
temperature  profiles  and  the  transmittance  model  to  perform  simulations.  First, 
linear  radiances  are  computed.  We  select  700  cm~^  as  a  reference  wave  number  and 
use  the  k  value  corresponding  to  the  generalized  weighting  function  of  Ch.  7 
determined  in  Subsection  5.2.  The  700  cm~^  Planck  intensity  profile  based  on  the 
atmospheric  (p  <  p,)  temperature  field  is  fitted  to  a  fifth-degree  polynomial. 
Then  from  Eq.  (2.17),  the  coefficients  of  the  polynomial  are  determined. 
Finally,  linear  radiances  are  computed  from  Eq.  (2.16). 

The  simulated  radiances  are  computed  from  Eq.  (2.1)  using  the  weighting 
functions  that  are  computed  from  the  transmittance  model.  Figures  9a- 1  show  the 
distribution  of  the  linear  radiances  vs.  simulated  radiances  for  the  12  channels 
used  in  the  temperature  retrieval.  We  see  that  t’-'e  correlation  for  each  channel 
is  significant,  except  Ch.  16.  Susskind  et  al .  (1983)  also  noted  that  the 


33 


WEIGHTING  FUNCTION 


34 


"atmospheric  response  of  Ch.  16  Is  very  diffuse."  This  may  be  due  to  the  fact 
that  the  weighting  function  of  Ch.  16  peaks  near  the  tropopause  at  about  200  mb 
(Fig.  8).  In  the  present  retrieval,  the  data  involving  Ch.  16  are  not  used.  The 
radiance  correlations  for  Chs.  6,  7,  and  13  are  particularly  significant, 
presumably  due  to  the  narrow  form  of  weighting  function  of  these  channels. 

We  then  perform  the  regression  analysis  for  each  channel  in  the  form 

-  a,  ♦b.R,.  ♦c.Rf,.  (5.1) 

where  R^]^  and  R^,  denote  linear  and  simulated  radiances.  Table  1  lists  the  values 
of  a^,  b^  and  c^  along  with  the  root -mean- square  (rms)  differences  between  the 
model -computed  and  fitted  values  of  R^^.  The  rms  differences  between  linear  and 
simulated  radiances  are  less  than  1.3  mW/mVsr/cm~^ ,  indicating  that  the  linear 
and  simulated  radiances  are  highly  correlated.  The  coefficients  a^,  bj^  and  c^ 
depend  on  the  functional  form  of  weighting  function,  the  reference  wavenumber 


Table  1  Coefficients  in  the  correlation  between  the  simulated  and  linear 
radiances  and  root-mean-square  differences  between  computed  and 
linear  radiances. 


Gh  • 

a 

b 

c 

rms  diff. 
(mW/mVsr/cm~^) 

1 

101.22 

-3.16 

0.0407 

0.62 

2 

-16.66 

1.35 

-0.0016 

0.37 

3 

-70.85 

3.46 

0.0214 

0.60 

4 

-9.19 

1.68 

-0.0110 

0,91 

5 

10.96 

0.75 

0.0004 

0.74 

6 

-2.53 

1.21 

-0.0016 

0.28 

7 

-3.51 

1.37 

-0.0024 

0.77 

13 

51.30 

51.63 

-12.02 

1.29 

14 

46.65 

51.26 

-14.92 

1.26 

15 

42.41 

45.68 

-33.45 

1.22 

17 

31.26 

303.30 

-827.9 

0.86 
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and  the  degree  of  polynomial  for  the  linear  radiances.  These  coefficients  are 
insensitive  to  temperature  profiles,  however.  The  set  of  coefficients  shown  here 
only  represents  an  example  of  the  correlation  between  simulated  and  linear 
radiances.  Note  that  different  reference  wave  number,  sharpness  Index,  and 
polynomial  will  produce  different  correlation  equations.  Further  studies  are 
required  to  search  for  the  optimal  correlation. 
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Section  6 


RETRIEVAL  RESULTS  AND  DISCUSSIONS 


Using  the  BUAN  collocated  satellite  and  radiosonde  data,  we  have  performed 
retrievals  to  investigate  the  applicability  of  the  DI  to  satellite  data.  The 
data  consists  of  1000  profiles  selected  for  this  study  which  were  collected 
between  April  3  and  April  7,  1988,  and  these  profiles  encompass  the  entire  globe. 
Based  on  the  N*  values  given  in  the  BUAN  data,  the  profiles  have  been  separated 
into  three  groups:  clear  (N*  -  0),  partly  cloudy  (0  <  N*  <  1)  ,  and  overcast  (N* 
>  1).  Table  2  lists  the  geographical  distribution  of  the  selected  data.  More 
than  half  of  the  cases  are  clear.  About  20X  of  the  cases  are  identified  as 
partly  cloudy.  About  33X  of  the  cases  occur  over  the  northern  midlatitude  area 
(60*N-30*N),  while  about  SOX  of  the  cases  occur  over  the  oceans. 

Temperatures  are  computed  for  each  peak  pressure  level  of  the  weighting 
functions  between  17  and  100  mb.  There  are  eleven  layers  between  these  levels. 
For  each  of  the  ten  layers  between  28  and  1000  mb,  two-point  averages  are  applied 
to  obtain  layer-mean  values.  The  top  layer  between  17  and  28  mb  has  not  been 


Table  2  Geographical  distribution  of  selected  data. 


Geographical 

Zone 

Number  of 

Profiles 

Clear 

Partly  Cloudy 

Overcast 

Total 

90*N-60*N 

168 

36 

36 

240 

60*N-30’N 

130 

91 

101 

322 

30*N-30*S 

173 

47 

56 

276 

30*S-60*S 

66 

10 

36 

112 

60*S-90*S 

23 

13 

14 

50 

ocean 

472 

158 

171 

801 

land 

88 

39 

72 

199 

Total 

560 

197 

243 

1000 
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selected  for  presentation  because  of  the  large  tentperature  differences  in  this 
layer.  The  large  temperature  differences  are  due  to  errors  in  polynomial  fitting 
for  p  <  17  mb.  The  purpose  of  averaging  is  to  simulate  the  vertical  resolution 
of  the  input  temperature  field  for  weather  prediction  models.  These  models  are 
usually  Initialized  using  layer-mean  rather  than  level  temperature  values.  In 
addition,  differences  between  retrieved  and  sounding  values  are  smaller  based  on 
layer-means  than  based  on  level  values  because  of  the  smoothing  effect  of  the 
averaging  process. 

Differences  between  retrieved  and  sounding  layer-mean  temperature  values 
are  computed.  The  sounding  values  have  been  obtained  by  interpolation  of  the 
BUAN  temperature  values  of  the  15  mandatory  and  25  significant  levels.  However, 
data  are  not  available  at  some  levels.  For  example,  Table  3  lists  the  number  of 
clear  cases  for  each  layer  and  the  geographical  zone  where  sounding  values  are 
available.  By  comparing  the  first  column  of  Table  2  with  Table  3,  it  is  noted 
that  about  lOX  of  the  cases  in  the  troposphere  are  data -void.  The  reason  for 
these  missing  data  is  unknown.  Computation  of  differences  can  only  be  carried 
out  at  those  levels  where  sounding  values  are  available. 

Simulation  retrievals  are  carried  out  using  the  linear  radiances  as  input. 
The  purpose  is  to  examine  the  effects  of  errors  in  polynomial  fitting  on  the 
retrieval  performance.  Linear  radiances  have  been  generated  according  to  the 
method  described  in  Subsection  5.4.  Theoretically,  if  the  Planck  Intensity 
profile  based  on  the  soundings  can  be  represented  by  a  polynomial  in  an  exact 
manner,  then  the  retrieved  temperature  values  will  be  exactly  Identical  to  the 
soundings.  Because  the  Planck  Intensity  profile  cannot  be  exactly  expressed  in 
terms  of  a  polynomial,  retrieved  values  differ  from  the  soundings.  Table  4  lists 
the  mean  differences  for  each  layer  and  the  geographical  zone.  All  clear  and 
partly-cloudy  cases  are  included.  For  each  geographical  zone,  vertical 
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Tabic  3  Nuiter  of  clear  casca  for  each  layer  and  seographical  lone  uhere  aoundfng  vatuea  are  available. 
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variation  of  mean  difference  with  pressure  shows  a  systematic  pattern.  In 
general,  positive  values  occur  near  the  surface.  Groups  of  negative  and  positive 
values  occur  alternately  in  the  lower  troposphere,  upper  troposphere,  and 
stratosphere.  Clearly  this  is  due  to  systematic  errors  in  the  polynomial 
fitting.  The  largest  positive  differences  (-6  K)  occur  between  375  and  575  mb 
over  polar  regions.  These  large  values  are  associated  with  the  location  of  the 
tropopause.  The  polynomial  fitting  is  unable  to  resolve  the  sharp  inversion  near 
the  tropopause.  Since  the  differences  are  due  purely  to  the  mathematical 
representation  of  the  Planck  intensity  profile  and  are  systematic,  we  may 
subtract  these  mean  differences  from  the  retrieved  valuss  to  reduce  errors. 
Further  studies  are  required  to  relate  the  mean  differences  to  the  properties  of 
Planck  intensity  or  radiance  profiles. 

Comparison  of  the  results  from  different  radiance  input  reveals  the  effects 
of  adjustment  errors  and  instrument  noises.  Simulated  radiances  have  been 
generated  using  the  BUAN  temperature  and  humidity  profiles  based  on  Eq.  (2.1), 
and  are  adjusted  using  Eq.  (5.1).  Similarly,  BUAN  satellite  data  are  first 
converted  from  brightness  temperature  to  radiances  and  then  adjusted  using  Eq. 
(5.1).  Figure  10  shows  the  vertical  profiles  of  layer-mean  rms  temperature 
differences  for  retrievals  using  linear,  simulated  and  BUAN  radiances  as  inputs. 
All  clear  and  partly  cloudy  cases  are  included.  The  same  mean  difference  values 
given  in  Table  4  have  been  subtracted  from  the  retrieved  results  according  to 
latitudinal  zones  for  all  three  input  sets.  The  rms  differences  for  all  three 
sets  are  about  2  K  between  100  and  700  mb.  The  larger  differences  (-2.2  -  3  K) 
for  p  <  100  mb  are  mainly  due  to  inversions  near  the  tropopause  for  tropical 
areas.  The  rms  differences  for  retrieval  using  linear  radiances  are  less  than 
2  K  for  p  >  100  mb.  These  differences  are  caused  by  the  residual  effects  of 
polynomial  fluctuation  around  the  sounding  profile,  after  mean  differences  have 
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Fig.  10  The  vertical  profiles  of  rms  differences  from  temperature  retrievals 
using  linear,  simulated  and  BUAN  radiances  as  input. 


47 


been  subtracted.  The  ms  differences  for  slaulated  radiances  are  slightly  larger 
(<  0.5  K)  than  those  for  linear  radiances.  The  Increase  In  rms  differences  Is 
due  to  uncertainties  in  the  correlation  denoted  In  Eq.  (5.1).  Finally,  the  ms 
differences  for  BUAN  satellite  data  are  also  slightly  larger  (<  1  K)  than  those 
for  simulated  radiances.  The  Increase  In  rms  difference  near  surface  results 
from  the  differences  between  the  BUAM  and  simulated  radiances  associated  with 
Instrument  noise  and  error  In  the  weighting  function  of  the  near-surface  channels 
(Chs.  6,  13,  and  7).  It  appears  that,  between  100  and  700  mb,  the  effects  of 
correlation,  instrument  noise  and  error  in  the  weighting  function  on  resulting 
ms  differences  are  smaller  than  the  effects  produced  by  polynomial  fitting. 

To  compare  the  perfomance  of  DI  with  the  BUAN  retrieval  scheme,  the 
vertical  profiles  of  ms  differences  for  the  two  schemes  are  shown  in  Fig.  11. 
The  BUAN  satellite  data  have  been  used  In  both  schemes.  Also  shovm  is  the  ms 
differences  for  the  BUAN  first-guess  values.  All  clear  and  partly  cloudy  cases 
are  included.  Except  for  p  >  800  mb,  the  ras  differences  for  the  BUAN  scheme  are 
between  1  K  and  2  K.  For  p  >  100  mb,  the  ms  differences  for  DI  are  larger  than 
the  differences  for  the  BUAN  scheme  by  less  than  1  K.  For  p  >  500  mb,  the  larger 
ms  differences  for  DI  are  caused  by  the  neglect  of  the  effects  of  humidity 
variation  near  the  surface  on  the  weighting  functions  of  Chs.  6,  7,  and  13.  For 
300  mb  <  p  <  500  mb,  accuracies  of  DI  and  of  BUAN  scheme  are  comparable.  For  p 
<  300  mb,  the  larger  ras  differences  for  DI  are  produced  by  uncertainties  in  the 
polynomial  expressions  of  temperature  values  near  the  tropopause.  To  improve  DI , 
we  must  incorporate  the  effects  of  humidity  in  the  retrieval  and  improve  the 
fitting  of  temperature  inversions  at  the  tropopause.  The  results  shown  in  the 
diagram  indicate  that  DI  is  potentially  capable  of  achieving  retrieval 
performances  comparable  to  the  BUAN  scheme.  The  DI  neither  requires  a  first- 
guess  profile  nor  uses  error  covariance  matrices  to  regulate  solutions.  The  ras 
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Fig.  11  The  vertical  profiles  of  rms  temperature  differences  from  the  DI  and 
BUAN  schemes,  using  the  same  BUAN  radiances. 
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differences  for  the  BUAN  first -guesses,  as  shown  In  Fig.  11,  range  between  2  K 
in  the  upper  troposphere  and  5  K  near  the  surface.  The  ms  differences  for  DI 
are  everywhere  less  than  those  for  the  BUAN  first-guesses.  Thus,  DI  can  be  used 
as  a  tool  to  obtain  the  first-guess  profiles  for  the  mlnlBuoi- variance  aethod. 

To  test  the  proposed  multi -spectral  cloud- removal  scheme  described  in 
Section  4,  we  have  also  carried  out  retrievals  using  DI  under  partly-cloudy 
conditions.  Figure  12  shows  the  vertical  profiles  of  ms  differences  for  partly- 
cloudy  retrievals  using  the  DI  and  BUAN  schemes.  Two  sets  of  radiances  are 
employed  as  inputs  in  DI .  One  set  is  the  BUAN  radiances  which  were  obtained  by 
the  BUAN  cloud- removal  method.  The  other  set  is  the  cloud- removed  clear 
radiances  from  the  multi -spectral  method.  The  ms  differences  for  BUAN  first 
guess  profiles  are  also  shown.  Ue  have  selected  about  1/3  of  the  partly-cloudy 
cases  in  which  the  differences  in  Ch.  7  radiances  for  the  selected  pair  of  spots 
for  cloud- removal  are  larger  than  20  mW/m*/sr/cm“^ .  For  the  remaining  partly- 
cloudy  cases,  the  asstusptlon  involving  the  same  cloud  properties  may  not  be 
satisfied  so  that  the  slopes  and  Gg^  will  contain  large  errors.  McMillin 
(1978)  also  noted  that  the  clear  radiances  detemlned  for  cases  with  larger 
differences  in  the  window  channel  radiance  are  expected  to  be  closer  to  true 
values.  Similar  to  Fig.  11,  the  ms  differences  for  BUAN  retrievals  are  between 
1  and  2  K  for  p  <  800  mb.  The  rms  differences  for  DI  using  BUAN  radiances  are 
close  to  those  for  the  BUAN  scheme.  In  fact,  between  300  and  500  mb,  the  DI 
yields  better  accuracy  than  the  BUAN  scheme.  The  larger  ms  differences  near  90 
mb  are  due  to  uncertainty  in  the  retrieved  tropopause  temperature  over  tropical 
areas.  By  comparing  ms  differences  for  DI  using  two  sets  of  radiances,  we  see 
that  the  present  multi -spectral  method  produces  only  slightly  larger  rms 
differences  than  the  BUAN  N*-method  (-0.5  K) . 
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Fig.  12  The  vertical  profiles  of  rms  temperature  differences  for  partly* 
cloudy  conditions  using  the  D1  and  BUAN  schemes.  Two  sets  of 
radiances  are  used  in  DI;  one  is  the  cloud-removed  radiances  by  the 
current  multi -spectral  cloud-removal  method,  while  the  other  is  the 
cloud-removed  radiances  by  the  BUAN  N*-method. 
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However,  Che  larger  ms  differences  from  Che  presenc  mechod  show  chac  Che 
presenC  mulci-speccral  cloud- removal  scheme  requires  furcher  improvemenCs .  One 
such  improvemenC  involves  Che  decerainaCion  of  che  slopes  .  Addicional 
infomaCion  is  required  Co  examine  how  close  are  che  cloud  propercies  of  che 
selecCed  pair  of  spoCs.  Ve  may  derive  Chis  informaCion  by  using  Che  collocaced 
HIRS  and  AVHRR  satellite  data.  BoCh  HIRS  and  AVHRR  channels  are  on  board  NOaA 
polar-orbiclng  saCellices.  AVHRR  has  five  imaging  channels  wich  very  high 
resolucion  (-1  km  for  HPRT  daca) .  Thus,  one  single  HIRS  pixel  conCains  abouC  500 
AVHRR  pixels.  From  Che  AVHRR  daCa,  che  sCaCisCics  of  cloud  CemperaCure  and  che 
radiaCive  propercies  may  be  deCemined  for  each  HIRS  pixel.  This  informaCion  may 
be  used  Co  deCemine  a  more  reliable  slope  chan  Che  adjacenc- pixel  meChod. 
The  idea  of  using  che  collocaCed  HIRS  and  AVHRR  daCa  has  been  proposed  by  Aoki 
(1980) .  This  is  an  area  requiring  furcher  research  efforc  in  conneccion  wiCh  che 
presenc  cloud  removal  scheme. 

BeCween  100  and  800  mb,  Che  ms  differences  for  DI  using  Che  mulci-specCral 
mechod  are  abouC  2  K.  This  perfomance  is  comparable  Co  many  ocher  schemes 
(e.g.,  Susskind  ec  al.  ,  1984;  Eyre,  1989).  Moreover,  Che  presenc  mulCi-specCral 
mechod  can  be  applied  Co  cases  involving  boch  land  and  ocean  surfaces.  The 
operaCional  N*-meChod  uses  Che  combined  4.3  ^m  and  microwave  channels  Co 
deCemine  N*  over  che  oceans.  Over  land,  only  4.3  /im  channels  can  be  used.  As 
a  comparison,  Table  5  llsCs  che  vercically  averaged  ms  differences  for  boch  DI 
and  BUAN  schemes  over  land,  ic  has  been  well  known  Chac  CemperaCure  reCrievals 
over  land  are  less  sacisfaccory  due  Co  Che  nonhomogeneiCy  of  surface  CemperaCure. 
The  cable  shows  ChaC  che  vercically  averaged  ms  differences  are  all  larger  Chan 
2  K.  Using  Che  same  BUAN  radiances,  Che  average  ras  difference  for  DI  is  larger 
Chan  chac  for  BUAN  scheme  by  0.37  K,  due  Co  che  humidlCy  effecCs  near  Che 
surface.  However,  using  DI ,  Che  vercically  averaged  ms  difference  using  Che 


52 


Table  5  Vertically  averaged  ms  temperature  differences  for  both  DI  and  BUAN 
partly  cloudy  retrievals  over  land. 


BUAN 

DI 

first-guess 

retrieved 

BUAN  rad. 

CH678  rad. 

4.95  K 

2.40  K 

2.77  K 

2.86  K 

multi -spectral  radiances  is  only  0.09  K  larger  than  that  using  the  BUAN 
radiances.  This  shows  that  the  performance  of  the  present  multi -spectral  cloud- 
removal  scheme  is  comparable  to  the  BUAN  N*-method  over  land. 
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Section  7 


COHCLUSIONS 


In  this  report,  we  have  presented  the  applicability  of  DI  to  temperature 
retrievals  using  measured  satellite  data.  Ve  first  review  the  fundamentals  of 
DI  and  formulate  the  radiative  transfer  equation  that  is  appropriate  for  an 
infinite  atmosphere.  The  convolution  theory  of  bilateral  Laplace  transform  is 
subsequently  used  to  derive  an  expression  for  the  local  Planck  intensity  value. 
The  Planck  intensity  in  real  space  can  be  exactly  expressed  as  a  linear  sum  of 
the  higher- order  derivatives  of  radiance  profiles.  The  inversion  coefficient  in 
each  term  of  the  summation  can  be  obtained  based  on  a  generalized  form  of 
weighting  functions.  The  higher-order  derivatives  of  radiance  profiles  can  be 
obtained  by  fitting  Che  radiances  for  a  finite  number  of  channels  using  a 
polynomial  function. 

The  satellite  data  set  used  in  the  present  study  contains  collocated 
brightness  temperature  and  radiosonde  data  which  have  been  collected  during  the 
period  of  the  BUAN  experiment.  The  brightness  temperature  data  include  both 
cloud- contaminated  and  clear  cases.  Cloud-removed  clear  brightness  temperatures 
obtained  from  the  operational  N* -method  are  also  included.  The  satellite  data 
set  also  contains  the  first-guess  and  retrieved  temperature  profiles. 

A  new  cloud  removal  method  has  been  developed  in  conjunction  with  DI.  This 
method  reconstructs  clear  radiances  from  cloud- contaminated  brightness 
temperatures  using  the  principle  of  the  N*-method.  The  present  method  uses  the 
radiances  of  HIRS  Chs .  6 ,  7 ,  and  8  to  determine  clear  radiances  and  surface 
temperature  simultaneously,  based  on  radiative  transfer  simulations.  In  this 
multi -spectral  method,  radiance  correlations  for  the  three  channels  are  first 
developed.  The  linear  relationships  between  partly  cloudy  radiances  for  the  two 
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sets  of  channels  are  subsequently  used  to  obtain  another  set  of  clear -radiance 
equations.  The  above  correlation  equations  are  solved  sinultaneously  using  an 
iterative  approach. 

We  have  selected  HIRS  Chs.  1-7  and  Chs.  13-17  to  carry  out  temperature 
retrievals  using  Dl.  The  large  radiosonde  data  set  along  with  a  transmittance 
model  are  employed  to  compute  a  set  of  weighting  function  profiles .  Each  profile 
is  then  fitted  to  the  generalized  form  to  obtain  the  sharpness  index  parameters 
and  peak  pressure  levels.  The  average  values  of  these  parameters  are  used  in  the 
retrievals.  The  strategy  in  determining  the  optimal  degree  of  polynomial  to  be 
used  in  retrievals  is  subsequently  discussed.  To  apply  DI,  which  is  a  linear 
Inversion  method,  to  a  nonlinear  temperature  retrieval  problem,  we  have  developed 
an  adjustment  scheme  to  derive  linear  radiances  from  satellite  data.  We  compute 
a  large  number  of  linear  and  simulated  radiances  based  on  forward  radiative 
transfer  calculations.  For  each  channel,  a  correlation  in  the  form  of  quadrature 
fitting  is  established  between  the  linear  and  simulated  radiances.  These 
correlations  are  used  to  convert  measured  radiances  to  simulated  radiances,  thus 
removing  major  nonlinear  components  from  the  data.  The  rms  differences  between 
the  computed  and  fitted  linear  radiances  are  generally  less  chan  1.3 
mW/mVsr/cm~^,  indicating  that  the  linear  and  simulated  radiances  are  highly 
correlated. 

Results  of  retrievals  and  their  physical  significance  are  presented  and 
discussed.  Comparisons  of  retrieval  results  using  linear,  simulated  and  BUAN 
radiances  illustrate  the  effects  of  radiance  differences  due  to  correlations, 
instrument  noise,  collocation  error,  and  model  error.  The  retrieval  results  from 
DI  are  comparable  to  those  from  the  BUAN  scheme,  which  has  rms  differences 
between  1  -  2  K.  Near  the  surface  and  tropopause,  rms  differences  using  DI  are 
somewhat  larger  because  the  humidity  effects  are  neglected  in  the  retrieval.  In 
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addition,  the  temperature  inversion  at  the  tropopause  cannot  be  properly 
represented  by  polynomials.  Further  studies  along  these  lines  are  needed.  The 
present  multi -spectral  cloud-removal  scheme  produces  retrieval  accuracies  that 
are  comparable  to  operational  schemes.  Further  improvements  of  the  cloud- removal 
scheme  will  require  the  utilization  of  collocated  H1R5  and  AVHRR  data.  Finally, 
we  demonstrate  that  the  multi -spectral  cloud-removal  method  can  be  particularly 
useful  for  the  determination  of  surface  temperature  over  land,  where  the 
application  of  microwave  techniques  may  be  questionable. 
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Appendix 

BUAN  RETRIEVAL  PROCEDURES 


Starting  with  the  9x7  array,  nunerous  clear  tests  are  perfomed  on  each  of 
the  63  spots  to  search  for  clear  spots.  Those  spots  which  pass  all  the  clear 
tests  are  defined  as  clear.  The  rest  of  the  spots  are  subject  to  cloud-removal 
process  using  the  N*  method  (McMillin,  1978) .  Certain  restrictions  are  applied 
in  the  usage  of  the  N*  method.  Those  spots  in  which  the  N*  method  is  applicable 
are  defined  as  partly  cloudy.  The  remaining  spots  are  defined  as  cloudy. 

The  63  spots  are  further  divided  into  six  miniboxes,  each  of  which  contains 
an  array  of  3  x  3  spots  (dashed  frames  in  Fig.  3.1).  A  minibox  has  a  footprint 
of  50  km  at  the  nadir.  For  each  minibox,  the  radiance  sets  for  the  nine  spots 
are  examined.  If  one  or  more  spots  have  been  declared  clear,  then  the  minibox 
is  defined  as  clear.  Only  one  set  of  clear  radiances  is  produced  for  each 
minlbox.  If  only  one  spot  is  found  to  be  clear  in  the  minibox,  the  radiances  of 
this  spot  are  selected  to  be  the  representative  clear  radiances  of  the  mlnibox. 
If  more  than  one  spot  in  the  minlbox  is  declared  to  be  clear,  these  spots  are 
subject  to  a  test  to  determine  the  spot  that  is  the  "clearest."  This  test 
Involves  comparison  between  predicted  and  observed  MSU  channel  radiances.  The 
predicted  HSU  channel  radiances  are  obtained  by  an  interchannel  regression 
equation.  The  spot  that  shows  the  least  difference  between  predicted  and 
observed  MSU  radiances  is  flagged  as  the  "clearest." 

If  only  partly  cloudy  spots  are  contained  in  the  mlnibox,  the  minibox  is 
defined  as  cloudy.  Again,  only  one  set  of  cloud- removed  clear  radiances  is  to 
be  produced  to  represent  the  mlnibox.  For  this  purpose,  the  radiances  of  all 
partly-cloudy  spots  are  subject  to  tests  to  determine  the  "best"  cloud- removed 
clear  radiances.  These  tests  include  the  Interchannel  regression  test  that  is 
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performed  to  determine  the  "clearest”  spot,  and  the  search  for  the  aost 
successful  N*  correction.  The  radiances  for  the  "best"  partly  cloudy  spot  are 
selected  for  the  partly -cloudy  minlbox. 

If  none  of  the  spots  in  the  minlbox  are  either  clear  or  partly  cloudy,  the 
minlbox  is  defined  as  cloudy.  In  this  case,  the  spot  with  the  highest  value  of 
Ch.  7  radiances  is  selected.  The  representative  "clear"  radiances  for  a  cloudy 
minlbox  are  the  average  of  all  radiances  of  the  nine  adjacent  spots. 

For  each  box,  there  are  six  sets  of  representative  clear  radiances, 
representing  a  tremendous  amount  of  data  volume.  Further  filtering  tests  are 
then  performed  to  reduce  the  data  volume.  Each  minlbox  is  flagged  to  be  "good," 
"redundant,"  "questionable"  or  "bad."  Thus  only  those  mlnlboxes  which  are 
flagged  "good"  are  used  In  the  subsequent  retrieval  analyses. 

To  retrieve  atmospheric  temperature  profiles,  the  TOYS  operational 
retrieval  algorithm  is  used.  This  program  uses  a  minimum- variance  simultaneous 
retrieval  technique  (Fleming  et  al.,  1986).  Theories  of  the  minimum  variance 
method  have  been  described  by  Houghton  et  al.  (1984).  The  main  part  of  this 
retrieval  algorithm  is  to  solve  the  following  matrix  equation; 

dx  -  CK^  (KCK*  +  E)-i  Ay,  (A.l) 

where  Ax  denotes  the  difference  between  the  final  solution  and  the  first- guess 
solution  vectors,  C  is  the  expected  error  covariance  matrix  for  the  first-guess 
vector,  K  is  the  kernel  matrix  (Sy^/^Xj),  of  which  the  elements  can  be  derived 
from  the  radiative  transfer  equation,  E  is  the  expected  error  covariance  matrix 
of  the  combined  measurement  and  forward  radiative  transfer  model  errors,  and  Ay 
is  the  difference  between  the  observed  and  model -generated  radiances. 

Determination  of  the  first-guess  profile  is  an  important  step  in  the 
minimum-variance  method  because  the  final  solution  strongly  depends  on  it. 
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Procedures  to  compute  the  first-guess  vector  are  based  on  a  sample  of  previously 
collocated  radiosonde  and  satellite  observations.  Twenty  sets  of  radiances  among 
a  28-day  data  set  collected  prior  to  the  observation  date  are  selected.  They  are 
closest  to  the  observed  radiance  vector  for  a  given  sounding.  The  selected 
radiances  are  averaged  to  derive  the  initial  estimates  of  first-guess  radiances 
for  the  retrieval.  The  twenty  associated  temperature  and  moisture  profiles 
derived  from  the  statistical  regression  equation  are  also  averaged  to  produce  the 
first-guess  solution  vectors.  The  28-day  data  set  Is  categorized  according  to 
latitude,  land/sea  and  day/night  flags,  and  is  updated  each  week. 
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